/*
 * Main.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include "Applications/Equilibrator.h"
#include "Assemble.h"
#include "Multiphase.h"
#include "ReactionSystem.h"
#include "EquilibriumSolver/EquilibriumSolver.h"
#include "KineticsSolver.h"

#include <iomanip>
#include "Applications/ReactionPath.h"
#include "Utils.h"
#include "Test/Activity/HkfModelTest.h"
#include "Test/TRISpH/TRISpH.h"

int main()
{
	TRISpH(25.0, 1.0);
	
	//TestHkfModel_WaterNaCl(25.0, 1.0);
}

//int main()
//{
//	double T = 80.0;
//	double P = 100.0;
//	
//	vector<string> knownSpecies = {"H2O", "Calcite", "Magnesite"};
//	
//	const vector<string> initialConditions = 
//	{
//		// Equilibrium-Initial Conditions
//		"Composition::H2O       = 55.508 moles",
//		// Pure-Initial Conditions
//		"Composition::Calcite   = 2.0 moles",
//		"Composition::Magnesite = 2.0 moles"
//	};
//	
////	vector<string> knownSpecies = {"H2O", "Na[+]", "Cl[-]", "Calcite", "Magnesite", "CO2(g)", "H2O(g)"};
////	
////	const vector<string> initialConditions = 
////	{
////		// Equilibrium-Initial Conditions
////		"ChargeBalanced", 
////		"Composition::H2O         = 55.508 moles",
////		"TotalComposition::Na[+]  = 1.0 moles",
////		"TotalComposition::Cl[-]  = 1.0 moles",
////		"TotalComposition::Ca[2+] = 1.0E-7 moles",
////		"TotalComposition::Mg[2+] = 1.0E-7 moles",
////		"Composition::CO2(g)      = 2.0 moles",
////		// Pure-Initial Conditions
////		"Composition::Calcite     = 2.0 moles",
////		"Composition::Magnesite   = 2.0 moles"
////	};
//	
////	vector<string> knownSpecies = {"H2O", "Calcite"};
////	
////	const vector<string> initialConditions = 
////	{
////		// Equilibrium-Initial Conditions
////		//"ChargeBalanced",
////		"Composition::H2O = 55.508 moles",
////		//"TotalComposition::Ca[2+] = 2.0 moles",
////		//"TotalComposition::HCO3[-] = 2.0 moles",
////		// Pure-Initial Conditions
////		"Composition::Calcite = 2.0 moles"
////	};
//	
//	vector<string> allSpecies         = AssembleSpecies(knownSpecies, {"aqueous"});
//	vector<string> equilibriumSpecies = RemoveSpecies(allSpecies, {"basis", "mineral"});
//	vector<string> kineticSpecies     = FilterSpecies(allSpecies, {"mineral"});
//
//	ReactionPath path(allSpecies, equilibriumSpecies, kineticSpecies, initialConditions);
//	
//	path.SetTemperaturePressure(T, P);
//	
//	double hour = 3600.0;
//	
//	path.Solve(0.0, 3 * hour, 1.0E-6, "data.dat");
//}

//int main()
//{
//	double T = 80.0;
//	double P = 100.0;
//	
//	vector<string> knownSpecies = {"H2O", "Calcite", "Magnesite"};
//	
//	const vector<string> equilibriumConditions = 
//	{
//		"ChargeBalanced", 
//		"Composition::H2O = 55.508 moles",
//		"TotalComposition::Ca[2+] = 2.0 moles",
//		"TotalComposition::Mg[2+] = 2.0 moles",
//		"TotalComposition::HCO3[-] = 4.0 moles",
//	};
//	
//	vector<string> allSpecies         = AssembleSpecies(knownSpecies, {"aqueous"});
//	vector<string> equilibriumSpecies = RemoveSpecies(allSpecies, {"basis"});
//	
//	Equilibrator equilibrator(allSpecies, equilibriumSpecies, equilibriumConditions);
//	
//	equilibrator.SetTemperaturePressure(T, P);
//	
//	equilibrator.Solve();
//	
//	equilibrator.PrintState();
//}


//int main()
//{
//	double T = 80.0;
//	double P = 1.0;
//	
//	vector<string> knownSpecies = {"H2O", "Na[+]", "Cl[-]", "CO2(aq)"};
//	
//	vector<string> allSpecies = AssembleSpecies(knownSpecies, {"aqueous"});
//	
//	vector<string> equilibriumSpecies = RemoveSpecies(allSpecies, {"basis"});
//	
//	Multiphase multiphase = AssemblePhaseSystem(allSpecies);
//	
//	ReactionSystem equilibriumReactions = AssembleReactionSystem(equilibriumSpecies); 
//	equilibriumReactions.SetIndexes(multiphase.GetSpeciesMap());
//	
//	EquilibriumSolver equilibriumSolver(multiphase, equilibriumReactions);
//	
//	equilibriumSolver.SetTotalCompositionMode();
//	
//	VectorXd uj(5); uj << 1.0, 55.508, 0.02, 0.0, 1.0; // Cl[-], H2O, HCO3[-], H[+], Na[+]
//	
//	equilibriumSolver.SetTotalComposition(uj);
//	equilibriumSolver.UseReducedFormulation();
//	
//	VectorXd n = VectorXd::Constant(allSpecies.size(), 1.0E-4);
//	
//	const SpeciesMap& map = multiphase.GetSpeciesMap();
//	
//	n[map["H2O"]]   = 55.508;
//	n[map["Na[+]"]] = 1.0;
//	n[map["Cl[-]"]] = 1.0;
//	
//	cout << multiphase << endl << endl;
//	cout << equilibriumReactions << endl << endl;
//	
//	cout << setprecision(6);
//	
//	auto diagnostic = equilibriumSolver.Solve(T, P, n);
//	
//	cout << "Iterations:" << diagnostic.iterations << endl;
//	
//	multiphase.PrintState(T, P, n);
//	
//	equilibriumReactions.PrintState(T, P, n, multiphase.Activities(T, P, n));
//}

//int main()
//{
//	double T = 80.0;
//	double P = 100.0;
//	
//	vector<string> knownSpecies = {"H2O", "Halite", "Calcite", "Magnesite", "CO2(g)", "H2O(g)"};
//	
//	const vector<string> equilibriumConditions = 
//	{
//		"ChargeBalanced", 
//		"Composition::H2O = 55.508 moles",
//		"Composition::Calcite = 2.0 moles",
//		"Composition::Magnesite = 2.0 moles",
//		"TotalComposition::Na[+] = 1.0 moles",
//		"TotalComposition::Cl[-] = 1.0 moles",
//		"Composition::CO2(g) = 1.0 moles"
//	};
//	
//	Equilibrator equilibrator(knownSpecies, equilibriumConditions);
//	
//	VectorXd n = equilibrator.Solve(T, P);
//	
//	equilibrator.PrintState(T, P, n);
//}

///// Equilibrator Test Case
//int main()
//{
//	double T = 80.0;
//	double P = 1.0;
//	
//	vector<string> knownSpecies = {"H2O", "Calcite"};
//	
//	const vector<string> equilibriumConditions = 
//	{
//		"ChargeBalanced", 
//		"Composition::H2O = 55.508 moles",
//		"TotalComposition::Ca[2+] = 2.0 moles",
//		"TotalComposition::HCO3[-] = 2.0 moles"
//	};
//	
//	Equilibrator equilibrator(knownSpecies, equilibriumConditions);
//	
//	VectorXd n = equilibrator.Solve(T, P);
//	
//	equilibrator.PrintState(T, P, n);
//}

/// ReactionPath Test Case
//int main()
//{
//	double T = 80.0; 
//	double P = 100.0;
//	
//	double t    = 0.0;
//	double dt   = 0.1;
//	double tEnd = 10.0;
//	
//	vector<string> knownSpecies = {"H2O", "Calcite"};
//	
//	const vector<string> conditions = 
//	{
//		"ChargeBalanced", 
//		"Composition::H2O = 55.508 moles",
//		"TotalComposition::Ca[2+] = 1.0E-7 moles",
//		"TotalComposition::HCO3[-] = 1.0E-7 moles",
//	};
//	
//	vector<string> allSpecies         = AssembleSpecies(knownSpecies, {"aqueous"});
//	vector<string> equilibriumSpecies = RemoveSpecies(allSpecies, {"basis", "mineral"});
//	vector<string> kineticSpecies     = FilterSpecies(allSpecies, {"mineral"});
//	
//	Multiphase multiphase = AssembleMultiphase(allSpecies);
//	
//	ReactionSystem equilibriumReactions = AssembleReactionSystem(equilibriumSpecies);
//	ReactionSystem kineticReactions     = AssembleReactionSystem(kineticSpecies);
//	
//	equilibriumReactions.SetIndexes(multiphase.GetSpeciesMap());
//	kineticReactions.SetIndexes(multiphase.GetSpeciesMap());
//
//	EquilibriumSolver eSolver(multiphase, equilibriumReactions);
//	
//	vector<EquilibriumCondition> eConditions = EquilibriumConditions(conditions, multiphase, equilibriumReactions);
//	
//	eSolver.SetEquilibriumConditions(eConditions);
//	
//	VectorXd n = InitialGuessEquilibriumComposition(conditions, T, P, multiphase);
//	
//	n[multiphase.GetSpeciesMap()["Calcite"]] = 2.0;
//	
//	eSolver.Solve(T, P, n);
//	
//	multiphase.PrintState(T, P, n);
//	
//	cout << endl << endl;
//	
//	equilibriumReactions.PrintState(T, P, n, multiphase.Activities(T, P, n));
//	
//	KineticsSolver kSolver(multiphase, equilibriumReactions, kineticReactions);
//	
//	ofstream file("data.dat");
//	
//	file << scientific << setprecision(6);
//	
//	file << setw(15) << left << "Time"; 
//	BOOST_FOREACH(const string& species, multiphase.GetSpecies()) file << setw(15) << species;
//	file << endl;
//	
//	while(t < tEnd)
//	{
//		cout << "Progress: " << 100.0 * t/tEnd << endl;
//		
//		kSolver.Iterate(t, tEnd, dt, T, P, n);
//		
//		file << setw(15) << left << t;
//		for(int i = 0; i < n.size(); ++i) file << setw(15) << left << n[i];
//		file << endl;
//	}
//}

/// Equilibrator Test Case
//int main()
//{
//	Equilibrator equilibrator({"H2O", "Na[+]", "Cl[-]", "CO2(g)"});
//	
//	const vector<string> conditions = 
//	{
//		"ChargeBalanced", 
//		"Composition::H2O = 55.508 moles",
//		"TotalComposition::Na[+] = 1.0 moles",
//		"TotalComposition::Cl[-] = 1.0 moles",
//		"Composition::CO2(g) = 2.0 moles"
//	};
//	
//	equilibrator.SetEquilibriumConditions(conditions);
//	equilibrator.SetTemperaturePressure(25.0, 1.0);
//	equilibrator.Solve();
//	equilibrator.PrintState();
//}





//int main()
//{
//	ReactionPath path({"H2O", "Calcite", "Na[+]", "Cl[-]"});
//	
//	double T = 25.0, P = 100.0;
//	
//	map<string,double> n;
//	
//	n["H2O"]     = 55.508;
//	n["Calcite"] = 2.0;
//	n["Na[+]"]   = 1.0;
//	n["Cl[-]"]   = 1.0;
//	
//	path.Solve(100.0, T, P, n, "res.out");
//	
//	return 0;
//}

//#include <iostream>
//#include <iomanip>
//using namespace std;
//
//#include "OdeSolver.h"
//
//class MyODEProblem : public OdeProblem
//{
//public:
//	virtual void Function(double t, const VectorXd& y, VectorXd& f)
//	{
//		f[0] = +y[1];
//		f[1] = -y[0];
//	}
//
//	virtual void Jacobian(double t, const VectorXd& y, MatrixXd& dfdy, VectorXd& dfdt) {}
//};
//
//int main()
//{
//	MyODEProblem problem;
//	
//	OdeSolver solver(2);
//	
//	solver.SetControlSolution(1.0E-6, 0);
//	
//	VectorXd y(2); y << 0.0, 1.0;
//	
//	double t = 0.0, tEnd = 4.0;
//	double dt = 0.1;
//	
//	cout << setprecision(6) << fixed;
//	
//	while(t < tEnd)
//	{
//		if(!solver.Iterate(problem, t, tEnd, dt, y))
//			cout << "Error on iteration." << endl;
//		
//		cout << "t = " << t << '\t' << "y = " << y[0] << ", " << y[1] << endl;
//	}
//	
//	return 0;
//}


//#include <algorithm>
//#include <iostream>
//#include <iomanip>
//using namespace std;
//
//#include "Assemble.h"
//#include "EquilibriumSolver.h"
//
///// The application below compares the two formulations for the equilibrium problem solution: general and beta-invertible. 
///// For this case, the beta matrix is invertible, so the formulation beta-invertible can be applied.
///// While the general formulation achieves the max number of iterations (100), the beta-invertible formulation solves in one (1) iteration
//int main()
//{
//	double T = 25.0;
//	double P = 100.0;
//	const vector<string> knownSpecies = {"H2O", "Na[+]", "Cl[-]", "CO2(aq)"};//, "CO2(g)", "H2O(g)"};
//	const vector<string> allSpecies = AssembleSpecies(knownSpecies, {"aqueous"});
//	const vector<string> eqSpecies  = RemoveSpecies(allSpecies, {"basis"});
//	
//	PhaseSystem multiphase = AssemblePhaseSystem(allSpecies);
//	
//	cout << multiphase << endl;
//	
//	const SpeciesMap& map = multiphase.GetSpeciesMap();
//	
//	ReactionSystem reactionSystem = AssembleReactionSystem(eqSpecies);
//	
//	reactionSystem.SetIndexes(map);
//	
//	cout << reactionSystem << endl;
//	
//	cout << setprecision(6);
//	
//	EquilibriumSolver eqSolver(multiphase, reactionSystem);
//	
//	eqSolver.SetTotalCompositionMode();
//	
//	VectorXd uj(5); uj << 1.0, 55.508, 0.01, 0.0, 1.0;
//	
//	eqSolver.SetTotalComposition(uj);
//	
//	VectorXd n = VectorXd::Constant(allSpecies.size(), 1.0E-4);
//	
//	n[map["H2O"]]   = 55.508;
//	n[map["Na[+]"]] = 1.0;
//	n[map["Cl[-]"]] = 1.0;
//	
//	cout << setprecision(6);
//	
//	eqSolver.SetFormulation("general");
//	auto diagnostic = eqSolver.Solve(T, P, n);
//	
//	cout << "Iterations(General): = " << diagnostic.iterations << endl;
//	
//	cout << "n = " << endl << n << endl << endl;
//	
//	VectorXd a = multiphase.Activities(T, P, n);
//	
//	cout << "a = " << endl << a << endl << endl;
//	
//	VectorXd Ke = reactionSystem.EquilibriumConstants(T, P);
//	VectorXd Qe = reactionSystem.ReactionQuotients(a);
//	
//	cout << "log(Ke/Qe) = " << endl << (Ke.array()/Qe.array()).log() << endl << endl;
//	
//	eqSolver.SetFormulation("reduced");
//	diagnostic = eqSolver.Solve(T, P, n);
//	
//	cout << "Iterations(Reduced): = " << diagnostic.iterations << endl;
//	
//	cout << "n = " << endl << n << endl << endl;
//	
//	a = multiphase.Activities(T, P, n);
//	
//	cout << "a = " << endl << a << endl << endl;
//	
//	Ke = reactionSystem.EquilibriumConstants(T, P);
//	Qe = reactionSystem.ReactionQuotients(a);
//	
//	cout << "log(Ke/Qe) = " << endl << (Ke.array()/Qe.array()).log() << endl << endl;
//	
//	return 0;
//}

